; 29 Simpson's rule
; I = dx/3 * [f_0 + 4f_1 + 2f_2 + ... + 2f_{N-2} + 4f_{N-1} + f_{N}]
; define the sum factory rewrite to iter version (ex 31)
(define (sum_fac term next)
    (define (sum a b res)
        (if (> a b)
            res
            (sum (next a) b (+ res (term a)))
        )
    )
    (lambda (a b) (sum a b 0))
)
; sanity check
; (define int_sum (sum_fac (lambda (x) x) (lambda (x) (+ x 1))))
; (int_sum 1 100)
; break the sum into 2 parts.
(define (simpson f a b dx)
    (* dx (/ 1 3)
    (- (* 2 
        (+ ((sum_fac f (lambda (x) (+ x dx)) ) a b)  
            ((sum_fac f (lambda (x) (+ x (* 2 dx)))) (+ a dx) (- b dx)) )) 
        (f a) (f b)
    )
    )
)
(simpson (lambda (x) x) 0 10 0.1)

; 31 define prod_fac
(define (prod_fac term next)
    (define (prod a b res)
        (if (> a b)
            res
            (prod (next a) b (* res (term a)))
        )
    )
    (lambda (a b) (prod a b 1))
)
; sanity check
; ((prod_fac (lambda (x) x) (lambda (x) (+ 1 x))) 1 5)
; wallis calculating pi/4
(* 4. ( (prod_fac (lambda (x) (* (/ (- x 1) x) (/ (+ x 1) x))) (lambda (x) (+ x 2))) 3 1000) )

; 32 define accumulate_fac to get further abstraction
(define (accumulate_fac combiner null_value term next)
    (define (iter a b res)
        (if (> a b)
            res
            (iter (next a) b (combiner (term a) res))
        )
    )
    (lambda (a b) (iter a b null_value))
)
; you can use this to define sum_fac=(accumulate_fac + 0)
; and prod_fac = (accumulate_fac * 1)

; 33 even further abstraction: filtered_accumulate
(define (filter_accum_fac combiner null_value term next filter)
    (define (iter a b res)
        (if (> a b)
            res
            (if (filter a)
                (iter (next a) b (combiner (term a) res))
                (iter (next a) b res)
            )
        )
    )
    (lambda (a b) (iter a b null_value))
)

; sum the square of primes from a to b
; first copy prime?
(define (prime? n)
    (define (iter cnt)
        (if (< n (* cnt cnt))
            #t
            (if (= 0 (remainder n cnt))
                #f
                (iter (+ 1 cnt))
            )
        )
    )
    (cond ((= 1 n) #f)
        ((= 2 n) #t)
        (else (iter 2))
    )
)
(define sum_square_prime 
    (filter_accum_fac + 0 (lambda (x) (* x x)) (lambda (x) (+ 1 x)) prime?)
)
; sanity check
(sum_square_prime 1 5)

; prod of all positive int < n, s.t. gcd(i,n)=1
; first copying gcd
(define (gcd a b)
    (if (= 0 b)
        a
        (gcd b (remainder a b))
    )
)
(define (prod_relaprime n)
    (
        (filter_accum_fac * 1
                        (lambda (x) x)
                        (lambda (x) (+ 1 x))
                        (lambda (x) (= 1 (gcd x n)))
        )
        1 (- n 1)
    )
)
(prod_relaprime 10)

; 34 
; (define (f g)
;     (g 2)
; )
; (f f) = (f 2) which will not have results for 2 is not applicable


; 36 print the results in fixed-point
(define (fixed_point f x0)
    (display "\nFix Point Iteration\n")
    (define (iter x tol)
        (display x)
        (newline)
        (let ((next (f x)))
            (if (< (abs (- x next)) tol)
                next
                (iter next tol)
            )
        )
    )
    (iter x0 0.01)
)
(fixed_point cos 1)

; 37, 38 continued fraction with k-term truncation
(define (continued_frac n d k)
    (define (iter i res)
        (cond ( (= i 0) res)
            (else (iter (- i 1) (/ (n i) (+ (d i) res))))
        )
    )
    (iter k 0)
)
; calculating \phi = 1.618...
(/ 1. (continued_frac (lambda (x) 1) (lambda (x) 1) 100))
; calculating e
(define (d i)
    (cond ((= 2 (remainder i 3))  (* 2 (/ (+ 1 i) 3)))
        (else 1)
    )
)
(+ 2. (continued_frac (lambda (x) 1) d 100))
(+ 2. (continued_frac (lambda (x) 1) d 1000))
